% this function is used to test the biharmonic problem using argyris element 

% build the mesh
%[V,T] = mesh_init(25,2);  % trimesh(T,V(:,1),V(:,2),'Color','blue');
[p,e,t] = initmesh('squareg','Hmax',0.2,'Hgrad',1.99);
V = (p(1:2,:)'+1)/2; 
T = t(1:3,:)';
[T,E,ET,TE,EV] = build_fem_mesh(V,T);

[taux1,taux2,taux3] = directional_bary(V,T,1,0);
[tauy1,tauy2,tauy3] = directional_bary(V,T,0,1);

caseNum = 2;

% build the fem space
d = 5;
[dof_map_c0,n_dof_c0] = build_dof_map_c0(T,TE,ET,d);
[dof_map_c1,n_dof_c1] = build_dof_map_c1(T,TE,ET,d);
[M,W] = build_c1_weight(V, T, E, TE, d);


% using c0 scheme
dof_map = dof_map_c1;
dim = n_dof_c1;

% calculate stiff matrix with quradurature
x21 = V(T(:,2),1) - V(T(:,1),1);
x31 = V(T(:,3),1) - V(T(:,1),1);
y21 = V(T(:,2),2) - V(T(:,1),2);
y31 = V(T(:,3),2) - V(T(:,1),2);
areas = 0.5*(x21.*y31 - x31.*y21);

[qw,qp] = quad_rule(20);
[cr,desc,asce,div,J,K] = auxillary_mat(d);

L1 = qp(:,1); L2 = qp(:,2); L3 = 1- L1 - L2;
phi = vdm23(d,L1,L2,L3);
Mat1 = vdm23(d-1,L1,L2,L3);  % low order vdm
Mat2 = vdm23(d-2,L1,L2,L3);  % low order vdm
Id = diag(ones((d+1)*(d+2)/2,1));


rhs = zeros(dim,1);
m = (d+1)*(d+2)/2;  mm = m*m*size(T,1); 
S = zeros(mm,1);
Indx1 = zeros(mm,1); Indx2 = zeros(mm,1); 
pos = 1; % it may be accerate by cancel the loop.
for k = 1:size(T,1)   
    % 1. compute local system
    mat_loc = zeros(m,m);
    rhs_loc = zeros(m,1);
    
    % first derivatives on master elements
    u = d*de_cast_step(Id,d,taux1(k),taux2(k),taux3(k),desc);  % direction derivetives
    v = d*de_cast_step(Id,d,tauy1(k),tauy2(k),tauy3(k),desc);
    phi_x = Mat1*u;  %get the deriv B-form value
    phi_y = Mat1*v;

    % second derivatives on master elements
    uu = (d-1)*de_cast_step(u,d-1,taux1(k),taux2(k),taux3(k),desc);
    uv = (d-1)*de_cast_step(u,d-1,tauy1(k),tauy2(k),tauy3(k),desc);
    vv = (d-1)*de_cast_step(v,d-1,tauy1(k),tauy2(k),tauy3(k),desc);
    phi_xx = Mat2*uu;  %get the deriv B-form value
    phi_xy = Mat2*uv;
    phi_yy = Mat2*vv;
 
    % the right hand side value
    qx = V(T(k,1),1)*qp(:,1) + V(T(k,2),1)*qp(:,2) + V(T(k,3),1)*(1-qp(:,1)-qp(:,2));
    qy = V(T(k,1),2)*qp(:,1) + V(T(k,2),2)*qp(:,2) + V(T(k,3),2)*(1-qp(:,1)-qp(:,2));    
    val = test_fun(qx, qy, caseNum);
    
    % calculate the local matrix
    for i = 1:m
%          for j = 1:m
%             mat_loc(i,j) = (phi(:,i).*phi(:,j))'*qw*areas(k);
%          end
         mat_loc(i,:) = (phi_xx(:,i).*qw*areas(k))'*phi_xx + ...
             (phi_yy(:,i).*qw*areas(k))'*phi_yy;

         rhs_loc(i) = (phi(:,i).*val)'*qw*areas(k);
    end
    
    % compress to C1 d.o.f.
    mat_loc = W(:,:,k)'*mat_loc*W(:,:,k);  
    rhs_loc = W(:,:,k)'*rhs_loc;

%   assemble to the global system
    [i,j,s] = find(mat_loc); L = length(i);
    Indx1(pos:(pos + L - 1)) = dof_map(i,k);
    Indx2(pos:(pos + L - 1)) = dof_map(j,k);
    S(pos:(pos + L - 1)) = s;
    pos = pos + L;    

    rhs(dof_map(:,k)) = rhs(dof_map(:,k)) + rhs_loc;
end

K = sparse(Indx1(1:(pos-1)),Indx2(1:(pos-1)),S(1:(pos-1)),dim,dim);

% imposing previous c1 constrains
% [H,row_idx] = smooth_C1(dof_map,V,T,TE,ET,d);
% epsilon = 1e-3;
% c = lagrange22(K,rhs,H,zeros(size(H,1),1),epsilon,50,1e-6);


% imposing the boundary condition



c = K\rhs;

c0 = zeros(n_dof_c0,1);
for k = 1:size(T,1)
    c0(dof_map_c0(:,k)) = W(:,:,k)*c(dof_map_c1(:,k));
end
% c0 = c;

% calculate the error pointwisely
[vtx,tri,ph] = heval(dof_map_c0,V,T,c0,d);
p = test_fun(vtx(:,1),vtx(:,2), caseNum);
trisurf(tri,vtx(:,1),vtx(:,2),ph-p);